Modeling halo mass functions in chameleon f{R) gravity 
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On cosmological scales, observations of the cluster abundance currently place the strongest con- 
straints on f[R) gravity. These constraints lie in the large-field limit, where the modifications of 
general relativity can correctly be modeled setting the Compton wavelength of the scalar field to 
its background value. These bounds are, however, at the verge of penetrating into a regime, where 
the modifications become nonlinearly suppressed due to the chameleon mechanism and cannot be 
described by this linearized approximation. For future constraints based on observations subjected 
to cluster abundance, it is therefore essential to consistently model the chameleon effect. We analyze 
descriptions of the halo mass function in chameleon f[R) gravity using a mass- and environment- 
dependent spherical collapse model in combination with excursion set theory and phenomenological 
fits to A''-body simulations in the ACDM and f{R) gravity scenarios. Our halo mass functions 
consistently incorporate the chameleon suppression and cosmological parameter dependencies, im- 
proving upon previous formalisms and providing an important extension to A'^-body simulations for 
the application in consistent tests of gravity with observables sensitive to the abundance of clusters. 
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I. INTRODUCTION 

In ,f{R) gravity, the Einstein- Hilbert action is sup- 
plemented with a free nonlinear function f{K) of the 
Ricci scalar R [l|l, which when designed appropriately 
can contribute to produce the observed late-time accel- 
erated expansion of our universe [2n3]- /(^) gravity is 
formally equivalent to a scalar-tensor theory where the 
additional degree of freedom is described by the scalaron 
field Jr = d//di? '^-{Jl and the kinetic coupling vanishes 
in Jordan frame. Here, we parametrize our models by 
the scalaron field evaluated at the present background, 
I //JO I- The /fl field is massive, and below its Compton 
wavelength, it enhances gravitational forces by a factor 
of 1/3, increasing the growth of structure. Due to the 
density dependence of the scalaron's mass, f{R) gravity 
models may incorporate the chameleon suppression [8|- 
llO| , returning gravitational forces to Newtonian relations 
in high-density regions and making them compatible with 
Solar-System tests [Ulil- 

The enhanced gravitational coupling at low curvature 
and below the Compton wavelength can be utilized to 
place constraints on the f{R) modification. The transi- 
tion required to interpolate between the low curvature of 
the large-scale structure and the high curvature of the 
galactic halo [Il| as well as the comparison of nearby 
distances inferred from Cepheids and tip of the red giant 
branch stars in a sample of unscreened dwarf galaxies |13| 
set the currently strongest bounds on the background 
field, |/;^o| < 1*1 - (10^^ - 10-5), i.e., the typical depth 
of cosmological potential wells. Independently, strong 
constraints can also be inferred from the cosmological 
structure only. In this large-scale regime, the currently 
strongest constraints on f{R) gravity models are inferred 
from the analysis of the abundance of clusters, yielding 
a constraint of |/™| < 10"^ [li-El- 

It is important to note that the cluster-scale con- 



straints have been derived by relying on a linearized ap- 
proach of the f{R) modifications, assuming a linear rela- 
tion between the curvature fluctuation 5R and the field 
fluctuation SJr that is correctly described by the back- 
ground Compton wavelength of the scalaron. This ap- 
proach breaks down when |/flo| ^ IQ-^, where cluster 
scales are affected by the chameleon suppression. It is 
therefore important for comparison to future measure- 
ments to describe the observable quantities encompassing 
the chameleon effect. While A^-body simulations provide 
a great laboratory for the study of the chameleon mecha- 
nism |17H22 | , semi-analytic models need to be developed 
based on these simulations, in order to allow for a full 
exploration of the cosmological parameter space in the 
model comparison to observations |23l - [28| . 

In this paper, we develop and compare different 
prescriptions for modeling the halo mass function 
in chameleon f{R) gravity based on the mass- and 
environment-dependent spherical collapse for chameleon 
theories [2y] applied to f{R) models in combination with 
excursion set theory and phenomenological fits to A^- 
body simulations. Our descriptions of the halo mass 
function incorporate the chameleon mechanism and cos- 
mological parameter dependencies and show good agree- 
ment with A^-body simulations and previous fitting for- 
mulae for f{R) gravity without the need of introducing 
new fitting parameters for the chameleon modification. 
Thus, they are well-suited for complementing A^-body 
simulations for the consistent comparison of f{R) grav- 
ity to observations that are sensitive to the abundance of 
clusters. 

The outline of the paper is as follows. In i}IIl we review 
/(i?) gravity with a particular focus on the Hu-Sawicki 
model [ll|. We discuss the linearized and suppressed 
regimes and a description of the transition between them 
by an estimation of the thin-shell thickness Q . In ijllll we 
examine the evolution and formation of structure in f{R) 



gravity, in specific, through the spherical collapse model 
for chameleon theories and extended excursion set theory 
with a conditional moving barrier. We further give here 
details about the iV-body simulations employed. mV\ is 
devoted to the modeling of the halo mass function in 
chameleon f{R) gravity and the comparison of the dif- 
ferent approaches based on the spherical collapse model, 
excursion set theory, and phcnomcnological fitting func- 
tions to iV-body simulations. We conclude in fjVl 



II. ,f{R) GRAVITY 



where the background has been subtracted, i.e., Sfn ~ 

fniR) - !r{R), SR= R- R, Spm = Pm- Pm- 

The cosmological background is assumed to be spa- 
tially homogeneous and isotropic, where we describe the 
scalar metric perturbations of its Friedmann-Lemaitre- 
Robertson- Walker metric by ^ = (5goo/(2goo) and $ = 
6gii/{2gii). The relation of ^ to the matter density and 
6R is given by the modified Poisson equation |ll| 



V^* 



2k\ 

— Spn 



IsRifn) 



(7) 



In f{R) gravity, the Einstein-Hilbert action is supple- 
mented with a free nonlinear function of the Ricci scalar 

1 



5 = — / d^x^^ [R + /(i?)] + S^ {^^; g,,) , (1) 



where. 



= Stt G, Sm is the matter action with mat- 



ter fields ■0m, and we have adopted natural units. We 
specialize here to metric f{R) gravity, where the connec- 
tion is of Levi-Civita type and the modified Einstein field 
equations are obtained as usual from the variation of the 
action Eq. ([T]) with respect to the metric 5^1/, 






g tj.1^ -"^ tNiyfn. = K Tf_i^. (2) 



The scalaron fn = df/dR is the additional scalar degree 
of freedom of the model, characterizing the modification 
of the gravitational force. 

We further concentrate on the functional form of f{R) 
proposed by Hu & Sawicki \V\\ . 



f{R) = -fh^ 



ci {R/f 



C2(i?/m2)"+l' 



(3) 



Here, m^ = k^ Pmo/3 and overbars refer to background 
quantities. We require the modification f{R) to satisfy 
Solar-System tests [ll[ through the chameleon mecha- 
nism [81-1 IQJI and moreover, yield a Hubble parameter that 
matches the ACDM expansion history. This constrains 
the free parameters of the model ci, C2, and n. At high 



curvatures, C2 R ^ m^ and Eq. ([3]) simplifies to 



n+l 



•^^ ' C2 n R"- 



(4) 



with fno = fR{Ro) and ^0 denoting the present back- 
ground curvature. Furthermore, from requiring ACDM 
to be recovered when |//jo| ^- 0, we obtain 



ci -2 
-m 

C2 



— m^ = 2k^ 



PA. 



(5) 



The scalar field equation follows from the variation of 
the action Eq. Q with respect to the scalaron and in the 
quasistatic approximation, for |/fl| <C 1, becomes 



V'SfR 



1 



[SR{h 



■Spn 



(6) 



A. Linearized and suppressed regimes 

For large values of the background field compared to 
the typical depth of gravitational potentials, l/^ol 3> 
1^1 ~ (10^^ — 10^^), we can linearize the field equations, 
Eqs. (|ni) and ([7]), using the approximation 



5R: 



dR 

Wr 



S/r ^ 3m 5fR, 



(8) 



R=R 



where m is the mass of the scalaron evaluated at the 
background and Ac = 2Trm~^ is its Compton wave- 
length. Within this linearized approximation and in 
Fourier space, the solution to Eqs. ^ and ([7|) becomes 



fc^*(k) = -:- i-^ 



(9) 

where k ~ |k| is the comoving wavcnumbcr. From 
Eq. ^, it can be seen that at scales k '^ ma, gravi- 
tational forces are enhanced by a factor of 1/3. 

In the opposite limit, where l/^ol ^ |^| ^ (10~^ — 
10"^), using Eqs. (g]), ©, and ^, in the high-density 
regions, where nSp^a ^ S\7^6fR, the scalar field be- 
comes 



Ir ~ Iro 



i?o 



K^(Pm +4pa) 



(10) 



Hence, for p^ ^ pc and n > —1, we get /jj ~ 0, a sup- 
pression of the modifications and a return to Newtonian 
gravity. More specifically, in this case, 5R ~ K^Spnj and 
the Fourier transform of Eq. ([7]) recovers the standard 
Poisson equation 



kH{k) = -—a'Sp^^ik). 



(11) 



The transition between the linearized and suppressed 
regimes, described by Eqs. (O and pTj) . respectively, for 
a top-hat overdensity may be approximated by an esti- 
mation of the thin-shell thickness as we shall discuss in 
the following. 



B. Transition between spherically-symmetric shells 
of constant density 



In the thin-shell regime, for r G [i?Oi^TH], the scalar 
field is ii, [2^ 



Khoury & Weltman |8| derived an estimation of the ra- 
dial profile of the scalar field (^ in a spherically-symmetric 
top-hat overdensity of radius i?TH with constant inner 
and outer matter density p^ and pout, respectively. On 
the inside and outside of i?TH , the solutions of the scalar- 
field, (/9out and (/3in, minimize the effective scalar- field po- 
tential Vcs{'^) defined by the scalar- field equation 



n'/' = Kff(<<5) 



(12) 



with the tilde denoting the Einstein frame. VeS consists 
of the scalar-field potential V{i^) and a contribution from 
the coupling of ip to the matter components. Ref. [8| find 
that the distance that is necessary for ^p to settle from 
LpoMt to iy9in is approximately given by 






K <^out - ^h- 



TH 



6/3 



* 



N 



(13) 



where /3 is defined by the transformation of the Jor- 
dan frame metric g^i, to the Einstein frame metric g^j^ 
through 



e-2^«^5M- 



9i,u = e -"^-^g^u- (14) 

The Newtonian potential at the surface of the sphere is 



*N = 



M 



Sir i?TH 6 



PinR 



in-K-THi 



with mass M = Air PmRj'v and hence, we obtain 



■•TH 



AR 

i?TH 



1 y^o 



^m 



PinR 



2 
■TH 



In f{R) gravity, P = -1/V& and g^^ 
Thus, for l/fll < 1, Eq. dTBl) becomes 

AR 



(1 



3 /j?,an — /i?„< 



-Rth 



K^Pin 



p2 



(15) 

(16) 
fR)9^lv■ 

(17) 



The inner and outer solution of the scalaron minimizing 
Vcs{ip) is equivalent to Eq. (fTO|). i.e.. 



/. 



i?, in/out 



l + A^ 



Pin/outfl 



/I Ha. 



n+1 



/« 



flO, 



(18) 



where An/out = Pm,in/out(a = 1)/Pm(a = 1). There- 
fore, the transition between spherically-symmetric shells 
of constant density in f{R) gravity can approximately be 
described by 



AR 

i?TH 




n+1 



(19) 



ip{r) ~ (/3i„ + 



K/3 



7' 



2^°'' 



(20) 



where AR ~ Rth — Ro- Hence, the extra force F for a 
unity test mass at Rth becomes 



■f3Wip\, 



^f^ 



2/3' 



.GM 

'r^~ 



V -Rth/ 



i?, 



TH 



^AR.JARVrARV 
Rth \ Rth ) \ Rth ) 



(21) 



Note that as Rth > -Ro and i?o > 0, we have ARJRth G 
[0,1], which implies F G [0,2/3^Fn] and, in specific, 
F £ [0,Fn/3] for f{R) gravity, where Fjq is the New- 
tonian force. Hence, for a top-hat overdensity, Eq. (PT|) 
yields an interpolation between the suppressed regime in 
Eq. ((TT|) and the 1/3 enhancement of the gravitational 
force in Eq. ©, which is C° for Ai?/i?TH ^ and C^ for 
AR/Rth -> 1. 

In previous studies [2a, [2^, only the first term in 
Eq. ([21]) has been considered through the approximation 



F: 



2/3Vin(3|^, 
Kth 



(22) 



This slightly underestimates the efficiency of the 
chameleon suppression. When studying the structure 
formation in chameleon f{R) gravity through the imple- 
mentation of the thin-shell approximation in the spheri- 
cal collapse model in the following, we shall use the full 
expression Eq. ((2T|). However, we also study the case of 
introducing a constant fudge factor a in Eq. 



F ~ 2/3^ min 3a 



AR 

Rth 



1 



(23) 



to modulate the efficiency of the chameleon suppression 
and account for corrections of approximations such as 
sphericity [301 ^^id a top-hat overdensity [2J| to realis- 
tic structure formation. Ref. [2^ found that a factor of 
a w 1/2 yields good agreement with the difference be- 
tween the lensing and dynamical mass of dark matter 
halos measured in iV-body simulations of f{R) gravity. 



III. STRUCTURE FORMATION 

In the following, we study the formation and evolution 
of structure in the cold dark matter scenarios of ACDM 
and f{R) gravity. We begin by reviewing the linear 
growth of structure for ACDM and the quasistatic regime 
of f{R) gravity in ^UTM where due to Eq. © for f{R) 
gravity, the linear growth becomes a function of scale in 
addition to its time dependence. In i jHIBI we describe 



the spherical cohapsc model for f{R) gravity, discussing 
its application to excursion set theory in mil CI We ex- 
amine the role of the environmental density in §111 D I and 
give details on the iV-body simulations employed in our 
study in flnH 



A. Linear growth of structure 

In ACDM, combining the linearly perturbed Einstein 
field equations with energy-momentum conservation, one 
obtains the ordinary second-order differential equation 
for the evolution of the matter overdensity Ani(a, fc) in 
total matter gauge 



A" 



-flni{a) 



A' 



-fin, (a) An 



0, (24) 



where here and throughout the paper, primes de- 
note derivatives with respect to In a and fiml^) = 
H^flmd^^/H^- We define the linear growth function 
D{a) as 



Dia) _ A^{a,k) 



D{a{) A„,(ai,fc) 



at an initial scale factor ai <C 1 in the matter-dominated 
regime and solve for D{a) with the corresponding initial 
conditions D{a\) = ai and D'{a\) = ai. Note that in 
this paper, D(cl) shall always refer to the linear growth 
function assuming a ACDM cosmology. 

In j(K) gravity, Eq. (PS|) is altered due to the modifica- 
tion of the Poisson equation, Eq. ((71), where for large |//j|, 
an additional modification of the relation of the lensing 
potential ($ — ^) to the matter density fluctuation con- 
tributes through the rescaling of its dependency on the 
matter density by (1 -I- Jr)^^ ■ These modifications can 
be included in Eq. (pi]) as 



a: 



2 - -fi,n(a) 



AL 



3 1 — .g(a, fc) 
2 1 + /k 



where as in Eq. ^ 

$ + * 



1 



$-* 



3 A;^ + m?d?' 



fini(a)Am ~ 0, 
(26) 

(27) 



and correctly describe the time- and scale-dependent lin- 
ear growth function Dj(R)(a-,k) defined as in Eq. ((25)) 
on quasistatic scales [l^. Note, however, that at near- 
horizon scales, for scalar-tensor models like /(-R) grav- 
ity, the matter fluctuation obtained from combining 
the linearly perturbed Einstein equations with energy- 
momentum conservation, in general, deviates from the 
matter fluctuation inferred from the quasistatic descrip- 
tion Eq. (|26p [31] . Since this modification is small in 
/(i?) gravity and additionally, here, we are interested in 
the high-curvature regime, we can safely neglect this con- 
tribution and furthermore the lensing modification as we 
restrict to models where \]r\ ^ 1. 



B. Spherical collapse 

We study the formation of clusters in /(i?) gravity 
using the spherical collapse model. We approximate 
the dark matter halo by a spherically-symmetric top-hat 
overdensity of initial radius i?TH with a constant mat- 
ter density pin and pout on the inside and outside. In 
order to incorporate the chameleon suppression in the 
spherical collapse calculation, we follow [2al and imple- 
ment the thin-shell thickness estimator for the chameleon 
transition by Q described in ^11 Bl in the case of f(K) 
gravity. We introduce ^(a) to denote the physical ra- 
dius of the overdensity at a, where i,(a\) = aii?TH- Note 
that the nonlinear evolution of this overdensity causes 
^(a) to deviate from this simple linear relation at a > Oi. 
This deviation shall be denoted by the dimensionless vari- 
able y = ^(a)/ai?TH- Conservation of mass enclosed 
in the overdensity implies p^a^R^^i ^ PmC"^ and hence, 
P = pm/pm = y~^- 

From Eq. (fT5|) , it follows that the thickness of the thin- 
shcU is 



(25) f^ 



4^ 



I f I „3n+4 
I IRQ I a 

fim(^O^TH)' 



r2/h 



l + 4gt 



n+1 



' l+4gt ^ 



ycnv 



n+1 



4?^«^ 



(28) 



where we use the notation yh and 2/onv to refer to the 
inner and outer overdensities, the halo and its local en- 
vironment, respectively. 

In order to describe the evolution of j/h, we model the 
effective modification to Newton's constant as 



Goff — 



l + Flf 



G, 



(29) 



where i^(A^/^) is given by the thin-shell approximation 
in ijIIBI With this modification, the equation of motion 
of the spherical shell is given by [2^, [32| 






2pA)--r(l + J^)'5pm, (30) 



which with pin = V\^ yields 

y'i^ 2 - ^an(a) y(, + ifi„.(a)(l-f-F) {y^^ - l) j/h = 0, 

(31) 
where dots denote cosmic time derivatives. Note that in 
Eq. (PU)) , we have used the contribution of the modified 
force through Fdpm rather than through Fp^, which 
was used in [2g. The final expression for the evolution 
of j/h in Eq. pip , however, is in agreement with Eq. (35) 
of ^. 

We assume that the environment follows a ACDM 
evolution, which in Eq. (|30p is obtained in the limit 



A^/^ ^- or cquivalcntly, F — )• 0. Thus, 



v" 

yenv 



2 - -f7,„(a) 



^/cnv + 2^«^(") (j^eni - l) 2/env = 0. 

(32) 

Eqs. ([5T|) and ([5^ form a system of coupled differential 
equations, which we solve by setting the initial conditions 
at ai <C 1 in the matter-dominated regime, 

-, Oh/onv.i ,„„x 

yh/onva — t , (ddj 



yh/c 



3h/c 



(34) 



We use the ACDM linear growth function D(a) from 
Eq. psp to extrapolate initial overdensities to present 
time, defining an effective linear overdensity 



3h/o 



/(x;Ci 



h/c 



^(ai) 



Jh/e 



(35) 



C. Excursion set theory 

Excursion sets correspond to regions where the matter 
density smoothed over this region exceeds a given thresh- 
old, defining the regions where virialized structures are 
expected to have formed [33l - [39| . The smoothed matter 
density perturbation field over a region of radius i? is 

5(x,i?) = Jw{\^-y\;R)S(y)d^y 



W{k; i?)5ke"'''d-^k 



(36) 



where W{\x — y\; R) is a window function and S{x) = 
Pm(x)/pm — 1 is the matter density perturbation. 
W{k;R) and Su are the corresponding Fourier trans- 
forms. For now, we shall only consider the initial density 
perturbation field and use (5(x) to refer to it. We as- 
sume (5(x) to be Gaussian. The density fluctuation field 
is characterized by its power spectrum P{k), for which 
the variance is 

S{R) = a^iR) = ((52(x; R)) = / M^(fc; R)P{k)Sk. 

(37) 
Hence, given the power spectrum, one can interchange 5* 
and R as measure of the scale of spherical perturbations. 
For a sharp window function W{k] R) in k-space, an in- 
cremental step in the smoothed initial overdensity field 
(5(x; R) in Eq. (pG]) is attributed to the extra higher-k 
modes. In this case, the wavenumbers arc uncorrelatcd 
such that the incremental steps satisfy the Markov prop- 
erty, i.e., steps depend on the current value only and 
are independent of previous values. The increment is a 
Gaussian field with zero mean and variance dS such that 
(5(x; S) can be described by a Brownian motion in S with 
Gaussian probability distribution 



P(S,S)dS 



V2^ 



*"/2^d<5. 



(38) 



For a scale-independent linear growth function, deter- 
mining the growth of both S and a/S", the linear density 
field remains Gaussian at all times. This holds partic- 
ularly for a ACDM universe and from now on, S{x; R) 
shall refer to the extrapolation of the smoothed initial 
matter density perturbation to z = via the linear den- 
sity growth function D{a) from Eq. (|25|) [see Eq. ([35]) ]. 
Note that since the linear growth function is simpler to 
calculate for ACDM than for f{R) gravity and especially 
due to its scale independence, we shall always use the 
linear ACDM growth function D{a) to do this extrapo- 
lation. Hence, for f{R) gravity, the extrapolated matter 
density field and associated quantities should be inter- 
preted as effective quantities only. 

In this spirit, a spherical region of initial radius R is 
considered to have collapsed to a virialized object today 
or live in a larger region which has collapsed earlier if 
(5(x; > R) > 5c, where 5c is the effective collapse density^ 
which may be determined from the initial matter over- 
density causing a singularity in Eq. pip and extrapolated 
to present time via D{a). In f{R) gravity, this critical 
density is dependent on the mass of the top-hat over- 
density and the local environment, (5c(x; M, Jenv) with 
M sa 47rpino^TH- We show the mass dependency of 5c 
for different |/flo| and 5cm for collapse today in Fig. [T] 
In general, in addition to the top-hat mass and the envi- 
ronment, 5c is dependent on r^m and the redshift of the 
collapse Zc- 

In a ACDM universe, 5c becomes independent of mass 
and environment and thus, for a given Q,^^ at a given 
Zc, defines a flat barrier (5^. In this case, the fraction 
of mass enclosed in virialized dark matter halos of mass 
M > 47r R^pni,i to the total mass is [43 



T{M, z) = 



1 



V2TiS 



e 



= -«V2S _ -{S-5^f/2S 



dS. 



(39) 
This corresponds to the fraction of Brownian motion 
trajectories which have crossed S^ at S. The fraction 
of mass enclosed in halos of masses corresponding to 
S{M) £ [S, S + dS] that collapse at z = Zc is given by 
the Prcss-Schechter expression [4l| 



4){S,Zc)dS 



1 



D{0) 5'c 



V2^D{zc) S 



1 D{0) {5'c 



'(40) 

where (f>{S) describes the distribution of Brownian mo- 
tion trajectories that flrst cross the barrier D{z ~ 
0)6^/D{z = Zc) at S. 

In f{R) gravity the barrier 6c is no longer flat and 
becomes dependent on S and the environment embed- 
ding the collapsing halo. In order to characterize the 
environment, we follow [26| and study a top-hat sphere 
with density perturbation (5env(x;x), evolving according 
to ACDM and specified by the choice of radius x, which 
embeds i5(x;i?). The crossing probability conditional on 
the Brownian motion trajectory passing 5cnv at S^ is 
4>[5c{S,5cnv), S\5cnv, S^]. This probability needs to be 
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FIG. 1: The collapse density Sc in chameleon f{R) gravity 
predicted by the spherical collapse model for different |/hoJ at 
Zc = 0. The chameleon effect is incorporated via Eq. (|21|l and 
the thin-shell thickness Eq. (|28[) . Note that the predictions 
for (5c in f{R) gravity return to the ACDM value 5^ ~ 1.676 
when the environmental density fluctuation approaches the 
value of the halo overdensity. 



computed numerically, for which we use a code develo ped 
in [2al based on the algorithm of [42|- We refer to [2g 
for more details on this computation. In the following, 
we shall clarify the characterization of the environment. 



D. Environment 

The effect of the environment on the first-crossing dis- 
tribution depends on the definition of the radius x- This 
can be done, for instance, following the fixed-scale envi- 
ronment approximation of [26J, fixing a Lagrangian (or 
initial comoving) radius ^, defining % = ^. We shall 
adopt the value for the Lagrangian radius used in |26| . 
^ = %h~^ Mpc such that S'^ = cr|. For comparison to iV- 
body simulations or observations that do not differenti- 
ate between structures formed in different environments, 
the conditional first-crossing distribution of the moving 
barrier (t)\5r.{S, Junv), '5'|(5(inv, Si\ computed with the algo- 
rithm of [42| needs to be integrated over all environments. 
In order to do so, in the following, we shall denote the 
distribution of (5env characterized through <^ as P^[5cm), 
corresponding to the probability that the Brownian mo- 
tion trajectory passes through ^cnv at S^ never having 



crossed the collapse density 5^ at S < S^. This is given 
bylH 



Pd^^ 



V2^ 



■■Q{5^ - <5cnv) 



e~^^ - e~ 



26V 



(41) 



where O is the Heaviside step function. The 

environment-averaged first-crossing distribution becomes 



(0(5))env= / Pc-'^['5c(^,'5env),5|<5cnv,55]dJenv. (42) 

J —oo 

Note that this reduces to the unconditional first- 
crossing of a constant barrier 5^ at S, 



{<PiS))c 



1 



-{S^f/2S 



V2^ S 
when Sc{S, Jonv) = S^, for which 

(l)[6c{S,dcm),S\5cnv,S(^] 



(43) 



K- 



/2^{S - 5^)3/2 



2(S-Sf) 



(44) 

However, in general, {4>{S))cnv must be computed numer- 
ically. 

A more accurate approach for defining the radius x 
is taken in [27|, where the size of environments is de- 
fined by the Eulerian (physical) radius C. We shall adopt 
the value used in [231, C = 5/i^^ Mpc. We refer to [23] 
and in particular Fig. 2 therein for a comparison of the 
Lagrangian and Eulerian definitions for the environment 
and implications for the structures inferred from that for 
chameleon theories. The Eulerian overdensity at time t 

is HH 



ANL(t) 



1- 



6!^ 



1-S" 



This defines the barrier for the environment 



<5fnv(^^onv) = ^c 



AL 



PmV 



-l/^c 



(45) 



(46) 



for which the first crossing thereof determines the 5{t) 
that a spherical region containing mass M^nv must have 
in order to evolve into an Eulerian volume V at t. For 
a power-law matter power spectrum P{k) with index n^, 
this becomes 



ASd 



C 



8/i-iMpc 



3/Si 



3/(3+n,)S^ 



where 



(47) 



1 r°° 

^^(Menv) = S{0 = 2^ J k^Pik)W^kOdk (48) 



with Lagrangian radius ^ such that Mem = 47r^^pm/3. 
The first-crossing probabihty of the moving barrier 
^em{Si) in [S^tS^ + dS'^], PciYv{S(}dS^, corresponds to 
the probability that an arbitrary point is located in an 
environment, which will have an Eulerian radius C, at Zc 
and for which 6{t) e [Sf,„{S^),6f^,{S^ + dS^)]. 

We use an approximation of the probability distribu- 
tion of the Eulerian environment (5env by |45| and also 
used in [23, lH, 



instead of virial masses for f{R) gravity in our predic- 
tions with the approximate relation 



Pd^e 



^(^/2 



V27r 
X exp 






1- 



5c 



2 (1 - (5env/<5c)' 

where /3 = {(/8)^/6c/<Tg , oj ~ (Sc7 with 

d In 5*^ Us + 3 

d In il/pnv 3 



j/2-l 



(49) 



(50) 



E. A'^-body simulations 

Dark matter A^-body simulations of f{R) gravity pro- 
vide a great laboratory for studying the chameleon mech- 
anism. Here, we use simulation results of [1^ [2^ 
for the comparison to the semi-analytic modeling de- 
scribed in i jIIIBI through ijIIIDl These simulations are 
performed using a particle mesh code solving the qua- 
sistatic relations Eqs. ([6]) and ([7]). They cover the New- 
tonian and chameleon scenarios for each field strength 
|/^o| = 10-6,10~^10"■* with n = 1 and cosmologi- 
cal parameter values fixed to match WMAP 3-year re- 
sults, flA = 0.76, n^ = 1 - n^, h = 0.73, ?is = 0.958, 
and the initial power in curvature fiuctuations Ag = 
(4.89 X 10-5)2 at fc = 0.05 Mpc'^ Each set of sim- 
ulations consists of 10 realizations with each box size, 
Lbox = 64/1-1 ^p^^ 128/1-1 Mpc, 256/1-1 Mpc, and a 
total particle number of iVp = 256"^ placed on 128'^ do- 
main grids. During the simulation, the domain grids are 
progressively refined in regions where the local densities 
are sufficiently large to reach a predefined threshold. This 
causes the grid structure to efficiently follow the density 
distribution so that the high density regions can be re- 
solved better. For the identification of halos within the 
simulation and their associated masses, a spherical over- 
density (SO) algorithm (cf. [43) is used. Hereby, the 
particles are placed on the grid by a cloud-in-cell inter- 
polation and counted within a growing sphere around the 
center of mass until the required overdensity is reached. 
The particle masses contained in the sphere then define 
the mass of the halo. This process is started at the high- 
est overdensity grid point and hierarchically continued to 
lower overdensity grid points until all halos are identified. 

Note that we use the virial overdensity Avir obtained 
for ACDM to identify halos even in f{R) gravity in order 
to make a fair comparison between the different models. 
We estimate the error of using ACDM virial masses Afvir 



A^vir,/(fl) 
Afvir 



^virj(fl) 

Z-ivir 



-1/3 



(51) 



which becomes exact for a halo density scaling as p ~ 
j,-9/4 This radial dependence for p can be motivated 
for the self-similar secondar y in fall and accretion in both 
ACDM and f{R) gravity [li, E^. For our choice of 
cosmological parameters, in the case i^ = 0, we have 
Avir = 390 and for F = 1/3, the virial overdensity be- 



■JiR) 



309 [32[. Hence, the error is approxi- 



comes 

mately 8% in case of the full modification F = 1/3. As 
this is a very simplified estimate for the mass ratio in 
Eq. (|5ip . the full computation requiring mass and en- 
vironmental dependence of F and the exact radial halo 
profile, and due to its relative smallness compared to the 
overall modification, we chose to ignore this effect when 
comparing models of the halo mass function to the A^- 
body results in the following. 

Finally, note that recently, it has been shown [4^ 
that for symmetron models [49[, differences may appear 
between the scalar field distributions produced in A''- 
body simulations when assuming the quasistatic limit 
and when accounting for time derivatives of the scalar 
field, respectively. In the f{R) gravity simulations used 
here, the scalar field sits at the bottom of the effective 
potential and never changes sign. This is different from 
the symmetron model considered in [4^ , where the sign 
of the scalar field can be different in different regions af- 
ter the symmetry breaking. Hence, we do not expect the 
same magnitude in the deviations of the simulation re- 
sults for f{R) gravity models and that the small-scale 
structure is correctly described by the quasistatic ap- 
proximation in Eqs. (jS]) and ([7]). Furthermore, for f{R) 
gravity, numerical self-consistency checks have been con- 
ducted [l3l , supporting the assumption of the smallness 
of the time derivatives. A more rigorous analysis of the 
applicability of the quasistatic approximation remains to 
be conducted in future work. 



IV. MODELING THE HALO MASS FUNCTION 

Effects from f{R) modifications of gravity on halo 
properties have been studied in, e.g., [la, [23l - [25l . [29l . [32l . 
[so . [Slf . The enhanced abundance of clusters caused by 
the modification was used in |14l4l6| in comparison to 
observations to place an upper bound on the scalaron 
background value of |/ii;o| ^ 10-^. However, given 
the expected constraints, these analyses have been car- 
ried out in the linearized regime of f{R) gravity, where 
the approximation Eq. (|S]) is valid. With future mea- 
surements, constraints will penetrate into the chameleon 
regime and it becomes important to consistently incor- 
porate the chameleon effect on the observables. 

Here, we focus on describing the halo mass function 
in f{R) gravity. Thereby, we use the spherical collapse 



model, excursion set formalism, and fitting formulae that 
have been calibrated to ACDM and f{R) gravity iV-body 
simulations. We restrict to f{R) models with exponent 
n = 1 corresponding to the choice of n in the A^-body 
simulations described in mil El Our relations can be used 
to explore halo mass functions in the cosmological pa- 
rameter space beyond the parameter values used in the 
A'^-body simulations and hence can be applied to consis- 
tently constrain \fRo\ in the chameleon regime. 



A. Excursion set theory 

Having determined the first-crossing distribution 
(/)(S', Scnv) hi mil CI and its environmental average 
{4>{S, (5onv))onv m ijIII Dl where (j){S, Senv) dS describes the 
fraction of mass enclosed in halos of masses correspond- 
ing to S{M) G [S, S + dS], the halo mass function can be 
computed from 



dn{M) 
dM 



dM 



M 



{(t){S,Scnv)) 



cnvy/cnv 



dS 



dM 



dM. (52) 



We show the relative difference of the halo mass function 
predicted by the excursion set approach for f{R) gravity 
outlined in OTCl and TOIDI with respect to ACDM in 
Fig.H 



B. Sheth-Tormen halo mass function 

Sheth & Tormen [521 introduced a modification of the 
Press-Schechtcr expression for the first-crossing distribu- 
tion as a function of the peak-threshold ;/ = 5c/ vS, 



V (f){v) ~ A\l —av^ 



(53) 



Here, ^ is a normalization parameter, i.e., / dv ^(y) — 1, 
and a = 0.707 and -p = 0.3. Eq. ([5^ is designed such that 
the halo mass function 



AM M^^ 'am 



(54) 



matches results from A^-body simulations and the mod- 
ification can be motivated by excursion set theory with 
a mo ving barrier such as caused through ellipsoidal col- 
lapse [sals!]. 

We compute the halo mass function defined by the 
first-crossing distribution in Eq. (|53p for the differ- 
ent /(-R) models assumed in the A^-body simulations 
in ^HIEI (5c is determined through the spherical col- 
lapse model in ijHIBI becoming dependent on mass and 
environment and entering Eq. (|53p through the peak- 
threshold V. We compare our predictions for /(i?) gravity 
to their counterpart from ACDM in Fig. [21 showing the 
relative enhancements of the halo mass function caused 
by the /(-R) modifications in different local environments 



(5onv and the environment-averaged case assuming the Eu- 
lerian distribution of Jenv given in Eq. P5|) . In Fig. [31 
we also show results from using a = 1/2 in Eq. ([25)1 to 
increase the efficiency of the chameleon suppression de- 
termined by the thin-shell expression Eq. ([19]). 



C. Nonlinear PPF formalism 

Li & Hu ^3 introduce a nonlinear parametrized post- 
Fricdmann (PPF) description to determine the halo mass 
function for chameleon /(i?) gravity. They phenomeno- 
logically interpolate between the linearized and sup- 
pressed regimes by introducing a chameleon PPF transi- 
tion in the variance as 



4pf(A'^) = 



5;g)(A/) 



(M/Mth)^ 51/cdm(^'0 



1 + (Af/A/th) 



(55) 



where Afth and /i are calibrated to fit simulation results. 
Assuming the same initial conditions for the /(i?) and 

1/2 

ACDM models, the variance S J.^-. is determined from 



Eq. (j37p with the linear power spectrum 



-Pf (i?) (a> fc) = 



D 



fiR) 



(a, fc) 



Dia) 



P{a,k), (56) 



where the linear growth functions are derived from solv- 
ing Eqs. ([Ml) and ^. The PPF peak threshold in [H 
is then given by 



s!' 



fpPF 



4pf(A0' 



(57) 



which they subsequently use in the Sheth-Tormen 
expression Eq. (|53p to approximate the halo mass 
function. We simultaneously fit (Mth,M), where 

A?th (lO^I/ijol) Mq/H = Afth, to the enhancements in 
the halo mass function obtained from the A"-body simu- 
lations fo_r \fm\ = 10"^, 10-^ 10"^ described in ^iHEl 
finding A/th ^ 2.172 x lO^^ ^nd ^i. ~ 1.415. We show the 
PPF fits for the enhancements in the halo mass function 
and the corresponding peak-thresholds in Fig. [31 

Based on our results for 5c from the spherical col- 
lapse model in ilHIBI applied to the Sheth-Tormen halo 
mass function in i jIVBI as an alternative determination 
of i^ppF, wc suggest to generalize and redefine the PPF 
peak-threshold as 



J^PPF 



5c{M,Scnv) 

S^/^{M) 



(58) 



where (•)env denotes the environmental average. Note 
that 5c is determined using the linear ACDM growth 
function to extrapolate the initial overdensity associated 
with the collapse and 5*^/^ is the variance obtained for 
ACDM. In this definition, the chameleon transition is 
incorporated within (5c(Af, (5cnv) through the estimation 
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FIG. 2: Relative enhancement of the halo mass function in chameleon f{R) gravity with respect to the prediction for ACDM. 
The environmental dependence is illustrated using the collapse density 5c from Fig. [l] computed with the spherical collapse 
model in i jlll Bl and applied to the Sheth-Tormen fit for ACDM simulations Eq. (|53[l {left-hand pandel). These predictions are 
averaged over the Eulerian environment defined in ijIIIDI (dashed line) and compared to the excursion set prediction (solid line) 
{right-hand panel). Note that the A''-body results at the low-mass end are contaminated by the inclusion of subhalos, which 
are not identified and removed in the SO halo-finder employed. 



of the thin-shell thickness Eq. (|19p . The advantage of 
this approach is that it is theoretically well-motivated, 
that i^ppF may be determined without calibration of fit- 
ting parameters to simulation results, and hence, that 
it encompasses dependencies on cosmological parame- 
ters and can easily be applied to other chameleon the- 
ories. We compare the different approaches for comput- 
ing i^ppF in Fig. [3l finding a good qualitative agreement 
between them and supporting the functional shape sug- 
gested in the phcnomenological PPF interpolation for- 
mula Eq. dnsi). 



V. CONCLUSION 

We have studied the spherical collapse of a top-hat 
overdensity in f{R) gravity, taking into account the 
chameleon suppression of modifications in high-density 
regions. The chameleon mechanism is approximated by 
an estimate of the thickness of a thin-shell interpolating 
the scalaron field between the constant spherical halo 
overdensity and the constant spherical environmental 
density. We implement this thickness estimation to ap- 
proximate the nonlinear evolution of the spherical over- 



density and the initial overdensity associated with the 
collapse. The collapse density obtained by this proce- 
dure is environment- and mass-dependent. 

We use excursion set theory to obtain the halo mass 
function predicted by /(i?) gravity and compare it to re- 
sults from A'^-body simulations. We further apply the 
peak threshold predicted by chameleon f{R) gravity to 
the Sheth-Tormen fitting function for the halo mass func- 
tions of ACDM A^-body simulations, to describe the en- 
hancement of the f{R) halo mass function relative to its 
ACDM counterpart. Thereby, halo mass functions are 
predicted for different environments, where we assume 
an Eulerian environment distribution to estimate an av- 
eraged result. Introducing a fudge factor in the thin-shell 
thickness to account for oversimplistic assumptions and 
approximations in the derivation of the chameleon barrier 
and to modulate the efficiency of the chameleon suppres- 
sion, we can improve the description of the enhancement 
at the high-mass end of the halo mass function observed 
in the simulations. This fudge factor is also preferred 
in the description of the difference between the lensing 
and dynamical mass of dark matter halos inferred from 
simulations [29| . 

Finally, we compare our results to a nonlinear PPF fit, 
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FIG. 3: Comparison of predictions for the relative deviation in the halo mass function {left-hand panel) and in the peak 
threshold {right-hand panel) for chameleon f{R) gravity with respect to the ACDM prediction, derived with the spherical 
collapse model in ijIIIBI Results are averaged over the Eulerian environment described in ijIIIDI where the halo mass functions 
are computed using the Sheth-Tormen expression Eq. (|53p . The agreement with the A''-body simulations at the high- mass end 
can be improved by introducing the fudge factor a = 1/2 in the thin-shell thickness used in the spherical collapse computation, 
increasing the efficiency of the chameleon suppression. We also compare our results to the phenomenological PPF fit of [23l |. 
providing a theoretical motivation for the functional form assumed therefor. Note that the A'^-body results at the low-mass end 
are contaminated by the inclusion of subhalos, which are not identified and removed in the SO halo-finder employed. 



which introduces a description of the chameleon mecha- 
nism by interpolating the variance of the matter fluctua- 
tions between the hnearized f{R) gravity regime and the 
fuUy-suppressed hmit, corresponding to ACDM. We find 
that the peak threshold predicted by our environment- 
and mass-dependent spherical coUapse computations, av- 
eraged over the Eulerian environment, is in agreement 
with the peak threshold of the nonlinear PPF descrip- 
tion, supporting the functional form suggested for this 
phenomenological fit. While the PPF interpolation pa- 
rameters have been fitted to iV-body simulations using 
particular cosmological parameters, however, our deriva- 
tion of the peak threshold in f{R) gravity may be applied 
free of fitting parameters and it furthermore incorporates 
cosmological parameter dependencies. Hence, our results 
can be used to extrapolate simulations beyond the set of 
simulated cosmological parameters for the use in parame- 
ter estimation analyses for inferring constraints on f{R) 
gravity employing observations sensitive to the cluster 



abundance. 
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